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SUMMARY OF OVERALL PROGRESS 

The research project investigates the failure modes of 
composite materials as a result low velocity impact by simula- 
ting the impact with a finite element analysis. An important 
facet of the project is the modeling of the impact of a solid 
onto cylindrical shells composed of composite materials. Most 
computer models have been developed for laminated composites 
in the form of plates and under uniform displacement or load. 
The model under development will simulate the delamination 
sustained when a composite material encounters impact from 
another rigid body. 

For the year, the investigators installed computer 
equipment, tested the computer network and developed a finite 
element model to compare results with known experimental data. 
The model simulated the impact of a steel rod onto a rotating 
shaft. The rod scribes the surface and creates a plowed re- 
gion which is subsequently measured optically for the surface 
area. The surface area was then compared with the results of 
the three dimensional finite element analysis. 

Pre-processing programs ( GMESH AND TANVEL ) were developed 
to generate node and element data for the input into the three 
dimensional, dynamic finite element analysis code (DYNA3D). 
The finite element mesh was configured with a fine mesh near 
the impact zone and a coarser mesh for the impacting rod and 
the regions surrounding the impact zone. For the computer 
simulations, five impacting loads (35, 50, 85, 185 and 285g) 
were used to determine the time history of the stresses, the 
scribed surface areas, and the amount of ridging. The proces- 
sing time of the computer codes amounted from 1-4 days. for 
example, the calculations for the 85g impact onto the rotating 
shaft approximated 45 hours of central processor unit's time. 
The calculated surface area were within 6-12 percent, relative 
error when compared to the actual scratch area. 


I . INTRODUCTION 

Composite materials have been increasingly used for struc- 
tural materials because of their high strength-to-weight 
ratio. These newly developed materials are by no means ex- 
cluded from damage resulting from low-velocity impact. In 
some cases of normal use, these materials may appear to be 
sound when struck by a foreign object, but interior structural 
damage may result from the separation between laminae. The 
delamination would not be immediately known and could later 
expand in use. 


Delamination has caused significant concern in the design 
and analysis of composite materials and structures. The 
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failure is complicated because it involves geometric and 
material discontinuities (i.e. interlaminar cracks and 
variation of ply properties in the transverse direction), as 
well as the coupled mode I, mode II, and mode III fracture in 
a layered anisotropic material system. For reference the 
three types of relative movements of the cracking surfaces 
have been shown in Fig. 1. Mode I represents the symmetric 
displacements with respect to the x-y and x-z planes. The 
sliding or shear mode. Mode II, has symmetry with respect to 
the x-y plane and distorted with respect to the x-z plane. 
The tearing mode, Mode III, has displacements that are 
distorted with respect to both x-y and x-z planes. 

A target-impactor test and including some models have been 
developed to study the effects of composite materials under 
different dynamic loading conditions. From results of this 
test, important data on the stress-strain-delamination rela- 
tionship can be used to make further studies on the behavior 
of multidirectional composites to low-velocity impact. 

In this section, a brief summary of the relevant litera- 
ture on composite materials is presented on failures resulting 
from impact and the failure modes of composite materials from 
a fracture mechanics viewpoint. The review is not an exten- 
sive discussion of the subject but will serve as an eye-opener 
to the complicated nature of failure modes of composite mate- 
rials. 


A. Impact Induced Failure 

The behavior of multidirectional composites to dynamic 
loads, as a function of laminae thickness, geometry, and load, 
is an important indicator of their damage resistance. In many 
structural applications impact behavior is considered one of 
the most important criteria for material selection [2]. The 
emphasis for the development of tougher graphite fiber rein- 
forced composites has brought an increase in composite impact 
testing. The data obtained from these tests can be used to 
provide valuable information about the impact failure mecha- 
nisms and the variables which can affect the impact resistance 
of graphite-fiber-reinforced polymer matrix composites [3]. 

The approach in studying the response of composite materi- 
als to low-velocity impact is shown in Fig. 2. The three 
major steps of the approach are: (1) determination of impac- 
tor-induced surface pressure and its distribution, (2) deter- 
mination of internal stresses in the composite target caused 
by the surface pressure, and (3) determination of failure 
modes in the target caused by the internal stresses [4]. 
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Shivakumar, Elber and Illg [5] have analyzed circular 
composite plates under static equivalent impact loads. Their 
analysis was based on the minimum potential energy method. 
The ply failure region and mode were calculated by using the 
Tsai-Wu [6] and the maximum stress criteria [7] given in 
Appendix A. A step-by-step incremental transverse displace- 
ment procedure was used to calculate plate load and ply 
stresses. Their analysis predicted the following: 

1. ) The first failure mode was splitting and initiated 

first in the bottom-most ply. 

2. ) The splitting failure thresholds were relatively low 

and tended to be lower for larger plates than for 
small plates. 

3. ) The splitting damage region in each ply was elongated 

in its fiber direction; the bottom ply had the 
largest damage region. 

4. ) Larger plates developed larger splitting damage re- 

gions than smaller plates subjected to the same 
transverse load. 

5. ) First fiber failure threshold (load or energy) are 

higher for large plates over a limited range of fiber 
radius-to-thickness ratio ( r/t > 10). 

Delamination in composite materials usually follows a 
geometrical pattern from the first ply to the last as reported 
by Sankar and Sun (8). In their investigation of graphite/ep- 
oxy beams, they reported that impact loads will readily create 
two types of cracks: (1) a single longitudinal crack on the 
back surface symmetric about the impact point and (2) two 
parallel longitudinal cracks on either side of the impact 
point. Also, a damage area on the front side (or impact side) 
was always on either side of the longitudinal axis of symme- 
try. On the back side there is usually two damage areas on 
either side of the transverse line of symmetry; whereas in the 
center, delamination occurs in all four quadrants. A represen- 
tative C-scan of the damage area is shown in Fig. 3. 

The plastic deformation at the delaminated crack tip can 
be accounted for in a number of ways by elasto-plastic frac- 
ture mechanics theories. Mahishi and Adams [9] made studies 
on the initiation and growth of delaminations in central and 
single-edge notched [+45/0] graphite/epoxy laminates. The 
symbolism for the laminate code has been explained in Appendix 
C. A three-dimensional elasto-plastic, generally orthotropic, 
finite-element method developed by Monib and Adams [10] was 
used and has been summarized in Appendix B. They assumed that 
the edge delamination involves only matrix-dominated fracture, 
based on elastic and brittle behavior, and that the crack 
surface is parallel to the ply interface. 
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Their criteria for delamination was based on the local 
elastic strain energy release rate in the presence of 
plasticity and they accounted for the coupled Mode I, II, and 
III delamination. Once their energy release rate was 
evaluated it was compared with experimentally measured 
fracture toughness values and when they were equal, the nodes 
were separated to simulate the crack growth. The 
experimentally measured [11] fracture toughness values for a 



The elastic strain energy release rate in the presence of 
plasticity for the three modes of fracture were given as: 


G - 1 

1 lAA 

< r zz/ K *z - 

(F z 

- ”L 4a ' 2 / R zzl 

(1) 

G II " 1 

11 TKk 

‘ F x/ R zx - 

(F x - 

T zx Aa)2//R zx 1 

(2) 

111 

tF y /K zy - 

(F y - 

^y aa)2 / R zyl 

(3) 


where 


Aa 

AA 


F x' F y' F z 

K zz' K zx' K zy 

*zz' K zx' R zy 

x y T y 

°zz ' zx ' T zy 


incremental crack growth, 

area of the incremental crack, 

force components, 

total stiffness coefficients, 

reduced stiffness coefficients, and 

yield stress values. 


The total elastic strain energy release rate is then 


+ G 


II 


+ G 


III 


( 4 ) 


Their analysis demonstrated that the proposed elasto-plastic 
fracture criterion is very effective in predicting stable 
delamination crack growth in composite laminates having non- 
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linear matrix-dominated properties. The result did predict a 
very large area of delamination because of the coarse grid 
used in their analysis. 

Joshi and Sun [12] experimentally studied crack patterns 
in conjunction with numerical results under impact induced 
fracture in a [0e/90,-/0e] laminated composite plate. They 
established the role of 5 transverse shear stress in proximal 
and middle layer crack initiation. Delamination zones were 
constructed using metallographic techniques of the damaged 
area. Their study is important in understanding the behavior 
of each laminae under dynamic loading by observing crack path 
versus shear stress maximum. 

Figure 4 is a typical fracture pattern of a transverse 
section approximately 0.5 mm from the impact point which was 
observed by Joshi and Sun [12]. Figure 5 shows details of the 
finite element mesh used to calculate the shear stresses shown 
in Figures 6 to 9. They concluded that a crack in the upper 
layer appeared to start at point B (in Fig. 4) and moved to- 
wards point C, from which delaminated crack propogates. A 
strong correlation can be seen (Figs. 6 to 9) between the 
crack path in the upper layer and the locus of the shear 
stress maximum. 


B. Failure Node Analysis 


One of the underlying principles of fracture mechanics is 
that unstable fracture occurs when the stress-intensity factor 
at the crack tip reaches a critical value. By knowing the 
critical value of K_, K , and K_ at failure for a given 
material of a particular thickness and at a specific tempera- 
ture and loading rate, the designer can determine flaw sizes 
that can be tolerated in structural members for a given design 
stress level. Conversely, for an existing crack that may be 
present in the structure the design stress level can be safely 
determined [13]. 


The interfacial cracks in multidirectional composite are 
introduced during fabrication or processing operations. These 
defects may be caused by incomplete wetting or trapped air 
bubbles between layers, or by debonding of two laminae as a 
result of high-stress concentrations at geometric or material 
discontinuities [14]. Kuo and Wang [15] defined the stress- 
intensity factors of the singular stress field around a crack 
tip as 
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where y is the imaginary part of the oscillatory stress 
singularity and 3 is a form factor. For the eight-node 
quadrilateral element used by Kuo and Wang [15] (Figs. 10 
and 11), the values for y and 3 were 1/2 and 1/4, respec- 
tively. The variables a and x are the crack length and 
element width, respectively. Very good agreement between 
their analysis and those of other references was reported as 
shown in Figs. 12 and 13. 

Other investigators [15-19] have also used the same 
methods as Kuo and Wang [15] for the analysis of interfacial 
cracks in composites. Lakshminarayana , Murthy and Srinath 
[16] applied their analysis to laminated anisotropic cylin- 
drical shells. Their numerical comparisons showed that the 
model will yield accurate stress intensity factors from a 
relatively coarse mesh. The mass and stiffness matrices for 
the laminated composites can be found in [16] or [17]. The 
finite element approach to analyze interlaminar stresses 
between composite laminae has proved to be quite accurate. 

Ross, Malvern, Sierakowski, and Takeda [1] analyzed 
[0 5 /90 5 /0 c] glass/epoxy plates subjected to central local- 
ized impact. Their analysis was performed using the finite- 
element computer program SAPIV [20]. This program incorpo- 
rated a three-dimensional elastic-orthotropic-material model 
with dynamic analysis capability that allows the calculation 
of transverse shear stresses directly. 

The numerical results based on an impact velocity of 
45.7 m/s are summarized in Figs. 14-17. The top inter- 
laminar plane is the one closest to the impact side and the 
bottom interlaminar plane is the one farthest from the im- 
pact side. For an elastic solution, the shear stress r 
attains or exceeds the estimated interlaminar shear strength 
of 17.24 MPa at almost every station along the plane paral- 
lel to the x-axis (fiber direction of the top and bottom 
laminas) in both interlaminar planes, as shown in Figs. 14 
and 15. The shear stress along the line at 90° to the top 
and bottom fiber direction exceeds the interlaminar shear 
strength only locally near the impact point, shown in Figs. 
16 and 17. Their results were based upon elastic material 
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models that may 
methodology to 
models into the 


prove useful in extending the analytical 
provide insight into incorporating failure 
evaluative procedures. 


C. Finite Element Analysis 

There are many different computer codes presently avail- 
able for calculating a dynamic response due to impact. 
Basically, there are two general classes of algorithms for 
dynamic problems: Lagrangian and Eulerian. Each has pecu- 
liar advantages and disadvantages for various problems. The 
primary advantage of Lagrangian codes is their ability to 
accurately track material boundaries and interfaces. They 
are ideally suited for constitutive models that require time 
histories of material behavior. Materials at an interface 
are designated as slave and master, and procedures are built 
into programs to permit contact, separation, and sliding, 
with or without friction, between master and slave surfaces. 
In contrast, Eulerian codes are ideally suited for treating 
large distortions. For problems in which large distortions 
predominate, an Eulerian description of material behavior is 
necessary. Early finite element three-dimensional computa- 
tions were impractical from the standpoint of both cost and 
computer-storage capability. Hence, the pioneering solu- 
tions to impact problems were often obtained with two-dimen- 
sional Lagrangian codes, such as HEMP [21] and TOODY [21], 
and two-dimensional Eulerian codes, such as HELP [21]. The 
advent of third- and fourth-generation computers such as the 
CDC 7600 and CRAY machines has made three-dimensional compu- 
tations feasible and attractive. Consequently, a number of 
computer programs capable of solving the fully three-dimen- 
sional equations of continuum physics were developed. 
Three-dimensional Lagrangian codes include EPIC-3 [22], 
HEMP3D [21], and DYNA3D. Some of the three-dimensional 
Eulerian codes are HULL, TRIOIL/TRIDORF , and METRIC [21]. 

In 1976, impacts of both a nickel cylinder and a nickel 
sphere onto an aluminum plate at 1500 m/s were performed by 
using a two-dimensional computer program, EPIC (Elastic- 
Plastic Impact Calculations) [23]. Johnson [23] reported 
that that a triangular element is better suited to represent 
severe distortions than a quadrilateral element because a 
triangle resists a node to cross the opposite side of the 
triangle during large distortion. For impact problems in- 
volving severe distortions, quadrilateral meshes become 
entangled and may take on shapes resulting in negative 
volumes, thus rendering computations useless. 

One year later, in 1977, a three-dimensional code, 
EPIC-3, was developed for high velocity impact problems. It 
was successfully applied to normal impact of a nickel trun- 
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cated cone onto an aluminum plate of variable thickness at 
1500 m/s and oblique impact of an aluminum rod onto a rigid 
surface at 1000 m/s. The formulation in EPIC-3 was derived 
for tetrahedral elements subjected to large strains and 
displacements (22]. In 1981, Hutchings et al. [24] assumed 
that a target is a rigid-plastic solid and a projectile is a 
rigid sphere. Their result provided a surprisingly success- 
ful explanation of the observed variation of crater volume 
and projectile energy loss with impact angle and velocity. 
For impact problems the impact force and duration are not 
easily determined. In 1985, two simple and improved models 
incorporating energy and mass balances were successfully 
developed to calculate impact force and duration during low- 
velocity impact of circular composite plates [25]. 

An explicit three-dimensional finite element code, 
(DYNA3D) developed by Lawrence Livermore National Labora- 
tory, California was first released in 1976. The code has 
evolved considerably since 1976 and originated from the 
early finite element codes by Wilson [26] and HONDO [27]. 
The interface treatment for handling contact, sliding, and 
impact has been the most significant contribution of DYNA3D . 
The approach has proven useful in many applications, such as 
"Bar Impact on Rigid Wall" [28], "Impact of Cylinder into 
Rail" [29], and "Three-Dimensional Finite Element Impact 
Analysis of A Nuclear Waste Truck Cask" [30]. 

The above examples cited deal with dynamic impact pro- 
blems using numerical methods of low-velocity impact or 
indentation of surfaces. The objective of this study is to 
combine finite element simulation with experimental com- 
posite material failures. In the past, these two disci- 
plines have usually been treated as two separate branches of 
study. Hereby, a numerical finite element code, DYNA3D , has 
been employed to simulate the impact system and provide 
predictive post-impact information for future laboratory 
experimentation as well as real-time impact behavior. A 
large number of engineering problems, such as the one pre- 
sented here, are too complex to be solved by closed-form 
analytical methods. For such problems, the finite element 
technique can be used to approximate the structural behav- 
ior during and after impact of solids onto composite 
materials . 

D . Summa ry 

The interlaminar fracture of composite materials has 
been extensively studied under low-velocity impact. From 
the results of impact studies, a generalized understanding 
of the complex behavior of multidirectional composites can 
be obtained. The study of these materials is complex and 
therefore requires some highly advanced method of obtaining 
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pertinent information such as: (1) What stresses cause frac- 
ture? and (2) How are these stresses related to the geome- 
trical and physical properties of each laminae? 

One method of obtaining information about this complex 
behavior is by a simple impact test. From this test, useful 
information concerning the transverse strain capability or 
stress capability of any material can be obtained. It has 
been shown that the impact behavior of composite materials 
is a function of laminae geometry and its properties. This 
will directly influence the stress distribution between each 
laminae under dynamic loading. 

From calculated interlaminar stresses, failure criteria 
can be used to initiate a crack. Once the crack has been 
initiated, or if there is an existing crack, then the strain 
energy release rate or the stress-intensity factors can be 
compared to experimentally determined values as criteria for 
crack propagation. With the aid of an algorithm and an 
efficient model the stresses between each laminae can be 
found quite rapidly. 


II. METHODOLOGY 

A schematic of the impact system assumed in the numerical 
simulation is shown in Fig. 18. The schematic was selected 
for the analysis because the principal investigators have 
experimental data to compare with the results of the finite 
element analysis. The shaft rotates with an tangential velo- 
city of 1000 rpm about the Z axis. A steel rod containing a 
diamond impacts the rotating shaft with a downward speed due 
to gravity loading. Although the shaft geometry is axisymme- 
tric, the loading due to impact between the rotating shaft and 
the diamond head is not. Therefore, the problem is not amen- 
able to a two-dimensional approach and must be solved with a 
three-dimensional model. An existing three-dimensional dyna- 
mic finite element structural analysis code, DYNA3D [6,7], was 
used to analyze the impact. From the analysis, a measure of 
anticipated surface area can be predicted and compared with 
that observed in the laboratory, to determine whether the 
computer model is correct. 

A description of the computational process, which is 
employed herein for the impact simulation depicted in the 
previous section is shown in Figure 19. Programs were run on 
a Digital Equipment Corporation VAX 11/780 and the acquired 
VAX Station II/GPX computers. Two pre-processors, GMESH and 
TANVEL, were developed to prepare the information in a form 
usable by the main program, DYNA3D . Results of the computa- 
tions typically contain physical quantities and material con- 
ditions throughout the physical domain as a function of time. 
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The resulting output files were generally too massive to read 
alphanumerically . Hence, TAURUS 131], a post-processor, pro- 
vides graphical displays of items of interest. Each component 
of these tandem codes is described in more detail in this 
section . 

An interactive program, named GMESH, was developed for 
the trial study to generate preliminary node and element data 
for input to the analysis code, because DYNA3D does not yet 
have significant mesh generation capabilities. The flowchart 
of GMESH, which is fully interactive, is shown in Figure 20. 
The program was designed for prismatic, circular shafts or 
cylinders. The longitudinal axis of the circular shaft was 
designated as the Z-direction. The number and relative loca- 
tion of nodes must be the same for each cross-section. On VAX 
computers running the VMS operating system the following com- 
mands were used to compile and link GMESH: 


$ ASSIGN NODELE.DAT FOR006 
$ ASSIGN TANVEL.DAT FOR007 
$ FORTRAN GMESH 
$ LINK GMESH 
$ RUN GMESH 


After executing the foregoing commands, a screen message 
prompts for input as shown in Figure 20. After finishing 
interactive responses to GMESH, two data files were created. 
One is NODELE.DAT, the node and element raw data for DYNA3D ; 
the other is TANVEL.DAT, the input data for tangential velo- 
city generation program TANVEL. The purpose of developing 
this tangential velocity generation routine will be stated 
below. 

In this impact simulation both the diamond head (impact- 
or) and the shaft (target) were in simultaneous motion: one 
linear, and the other rotational. Because the current version 
of DYNA3D does not allow direct input of the angular velocity 
of the rotating shaft, Hallquist [32] suggested that the angu- 
lar velocity of the shaft be provided to DYNA3D by calculating 
the rotational tangential velocities of every nodal point on 
the shaft. A tangential velocity generation routine, TANVEL, 
was developed for this purpose. To execute TANVEL, the fol- 
lowing commands were entered: 


$ ASSIGN TANVEL.DAT FOR005 
$ ASSIGN TANVEL. OUT FOR006 
$ FORTRAN TANVEL 
$ LINK TANVEL 
$ RUN TANVEL 
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The commands then created an output file, TANVEL . OUT , contain- 
ing the tangential velocities of every nodal point. 

By concatenating NODELE.DAT and TANVEL. OUT into a single 
data file, DYNA3D.DAT, renders the input data for DYNA3D near- 
ly ready. Although NODELE.DAT and TANVEL. OUT are in the form 
required for DYNA3D input instructions, certain items specific 
to the analysis were added. These include material property 
data, termination time for computations, initial velocity of 
the diamond head, and the sliding interface definitions. 

The finite-element method models complicated problems in 
structural engineering and solid mechanics, with the resulting 
equations involving many thousands of degrees of freedom. 
Basically, the transient algorithms currently in use may be 
classified into two general classes: implicit and explicit 
[33]. Implicit algorithms tend to be numerically stable, 
permitting large time steps. However, the cost per step is 
high and storage requirements tend to increase dramatically 
with the size of the mesh. Explicit algorithms tend to be 
inexpensive per step and require less storage than implicit 
algorithms, but numerical stability requires that small steps 
be employed [34,35]. 

DYNA3D employs explicit algorithms to develop a three- 
dimensional finite element code whose equations are integrated 
by the central difference method. Hence, the time integration 
is only conditionally stable with respect to step size. Dur- 
ing execution, DYNA3D continuously monitors the step size and 
adjusts it to keep the calculation stable. Fifteen material 
models are available in DYNA3D which cover a wide range of 
material behavior, such as elastic, kinematic/isotropic 
elastic-plastic, soil and crushable foam, thermo-elastic-plas- 
tic, rubber, temperature dependent elastic-plastic hydro- 
dynamic, and so on [28]. 

According to the first few trial executions of DYNA3D , 
the rotating shaft, impacted by the diamond head, undergoes 
elastic-plastic deformations. Hence, material type 3, an 
elastic-plastic behavior with kinematic/isotropic work harden- 
ing, is exclusively used in this analysis. There are three 
types of sliding interfaces available in DYNA3D: sliding only, 
sliding-with-gaps , and tied. Because the diamond head re- 
bounds from the shaft at the end of the impact, a sliding- 
with-gaps interface, in which the contact-impact algorithm 
employed by the code allows gaps to open or close, is used 
between the diamond head and rotating shaft [36]. 

For implementation of DYNA3D , one surface of the inter- 
face, is identified as a master surface (i.e. diamond head), 
which was more coarsely zoned, and the other was identified as 
a slave surface (i.e. rotating shaft). Each surface is 
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defined by a set of three or four node quadrilateral segments, 
called master and slave segments, on which the nodes of the 
slave and master surfaces, respectively, must slide. Time 
zero was defined as the instant of impact of the diamond head 
and the shaft. 

The X-Y plane was the plane of symmetry with respect to 
the loading, so only one-half of the geometry in the positive 
2-direction was actually modeled, as shown in Figs. 21 and 22. 
The diamond head was initially specified to have the velocity 
resulting from a gravity loading, which is 2.53 m/s in the 
negative Y-direction. The angular velocity of the shaft is 
input in terms of the rotational tangential velocity X- and Y- 
components at every nodal point, as previously described. 

The command for executing DYNA3D was 


$ RUN DYNA3D 

DYNA3D will prompt: 

PLEASE DEFINE INPUT FILE NAMES OR CHANGE DEFAULTS: 

Then, the following command 
> DYNA3D I-DYNA3D 

caused DYNA3D to execute, where I-DYNA3D means that the input 
filename was DYNA3D.DAT. 

TAURUS, an interactive post-processor for DYNA3D , reads 
the output binary plotfiles generated by DYNA3D and plots 
contours, time histories and deformed shapes. TAURUS has 
three phases 

Phase 0: initialization. 

Phase 1: geometry display with contouring, 

Phase 2: time history processing. 

According to the termination time and time interval be- 
tween complete state dumps in the input control card, DYNA3D 
generated several binary plotfiles under the names of 
D3PLOT.DAT, D3PLOT01.DAT, D3PLOT02.DAT, etc. On VAX 

computers, the command 

$ RUN TAURUS 

causes TAURUS to prompt: 

PLEASE DEFINE INPUT FILE NAMES OR CHANGE DEFAULTS: 
After this, the command 
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> TAURUS G-D3PL0T 

allowed TAURUS to access all of the binary plotfiles and 
acted interactively, where G-D3PL0T defined that the binary 
plot filenames are D3PL0T.DAT, D3PLOT01.DAT, D3PLOT02.DAT, 
etc. The commands needed to effectively plot data are listed 
in Table 1. 


Table 1: Phase 1 and Phase 2 Commands in TAURUS 


Phase 1 Commands 

Phase 2 Commands 

CONTOUR 

ASCL 

DAM 

ASET 

DSF 

ELEMENTS 

M 

ETIME 

PHS2 

GATHER 

RAXY 

GTIME 

RAYZ 

MATLS 

RAZX 

MTIME 

RX 

NODES 

RY 

NRTIME 

RZ 

NTIME 

STATE 

OSCL 

TIME 

OSET 

TRIAD 

PHSl 

VIEW 

PRINT 

XSCALE 


YSCALE 


ZSCALE 



III. RESULTS AND DISCUSSION 

The finite element mesh used in this simulation contains 
1040 elements (1030 for the rotating shaft and 10 for the 
vertical steel rod containing the diamond tip, respectively) 
and 1339 nodes (1298 for the rotating shaft and 41 for the 
diamond head and attached vertical steel rod, respectively). 
Eight-node isoparametric solid elements were used for spatial 
discretization of the finite element mesh used in this simu- 
lation as shown in Figures 21 through 24. Figure 21 depicts a 
front wire-frame view of the impact system in the X-Y plane 
where the Z-axis goes directly into the paper and the origin 
is located at the center of the shaft. A side view is shown 
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in Figure 22. Instead of showing only one-half of the geome- 
try of the system, Figs. 23 and 24 are two different angle 
views of the entire system. 

Two primitive meshes in early simulations are shown in 
Figs. 24 and 25. In Figure 24, the entire mesh, containing 
840 elements and 1155 nodes, is symmetrical about the longitu- 
dinal axis. The configuration was inefficient and costly in 
numerical computation, because deformations of the shaft from 
impact were not considered near the top surface. A finer mesh 
was developed for the vicinity of the diamond head, where the 
effect of impact is predominant, includes 870 elements and 
1162 nodes as shown in Fig. 26. To reduce the amount of 
effort required to generate meshes, the interface treatment in 
DYNA3D permits sudden transitions in zoning with tied inter- 
faces. The more coarsely meshed side of the interface is 
recommended as the master surface, and the finer meshed side 
as the slave surface. Each master node should coincide with a 
slave node to ensure complete displacement compatibility along 
the interface. As indicated in Figs. 25 and 26 these slave 
nodes do not coincide with master nodes between the interface. 
Hence, incompatible displacements exist along these inter- 
faces. The total number of elements in Figs. 25 and 26 is 
much less than that of Fig. 21, which could save substantial 
computational cost. However, both meshes of Figs. 25 and 26 
were not employed in the simulation because they would result 
in incompatible displacement of nodes as indicated in the 
figures . 

An angular velocity of 1000 rpm (i.e. 105 rad/s in the 
negative Z-direction) is specified for the rotating shaft. 
The angular velocity was used for the tangential velocities at 
every nodal point in the shaft at time zero. The tangential 
velocity of each node obviously changes once the impact has 
started, but DYNA3D makes no provision for altering this velo- 
city with time. Fortunately, the duration of impact was only 
0.1 to 0.2 milliseconds, and it is assumed that the continuous 
change in velocities will not cause disastrous results for 
such a short response time. For the simulation, five loads of 
35, 50, 85, 185, and 285 grams were applied to the diamond 
head for five executions. 


A time history plot of mesh geometry was made to monitor 
the calculations and to determine that the diamond head has 
rebounded off of the shaft. For example, for an 85 gram 
applied load 105 microseconds elapsed before the diamond tip 
bounced from the surface. To ensure the impact event is com- 
pletely finished, the total numerical simulation time was 
extended to about 200 micro-seconds, which requires approxi- 
mately 45 hours of central processor unit time on a VAX 11/780 
super-minicomputer. This results in a time step being used 
for this simulation of approximately 7.0x10 seconds. A 
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small time step was used because the minimum element dimen- 
sion, which determines the critical time step size, was only 
0.39 millimeters. The finer the mesh, the smaller the 
required time step size. 

A sequence of deformed configurations is shown in Figures 
27 through 32. As time progresses, indentation proceeds fol- 
lowed by rebounding as shown in Figs. 27 through 29 and Figs. 
30 through 32, respectively. The vertical rod also undergoes 
some bending moment caused by frictional contact and shown in 
Figs. 29 through 32. Figure 33 shows another view of the mesh 
geometry at the end of impact. By using a nodal displacement 
magnification factor of ten, Fig. 34 provides a top view of 
the mesh geometry, which indicates an indented crater. A 
close-up view of the impact region of interest is shown in 
Fig. 35, where nodes and elements are numbered. At time zero 
node 101 on the shaft is directly impacted by node 1299, which 
is the diamond's lowermost point. 

Time history plots of displacement, velocity, and accel- 
eration in Y-direction (due to an 85-gram applied load) for 
nodes 1299 and 101 are shown in Figs. 36 through 38, respec- 
tively. Because the diamond and shaft interface allows slip- 
page, nodes 1299 and 101 were not always in contact with each 
other. At the end of impact the Y-displacement curve tangent 
of node 101 approaches the horizon as a consequence of the 
plastic deformation. From Fig. 37, the terminal vertical 
component speeds of nodes 1299 and 101 are 0.9 and approxi- 
mately 0.0 m/s, respectively. The nodes' accelerations oscil- 
late near zero (i.e. no contact force) at the end of the 
response (Fig. 38), which means that the diamond head and 
rotating shaft have completely separated from each other. 

The Y-displacement-time history for nodes 100, 101 and 
102 representing the left node, impact node, and right node, 
are respectively shown in Fig. 39. The calculations show that 
the indentor causes the right node to form a ridge as expected 
as well as the left node. 

The stress developed in the elements (841 and 851) sus- 
taining impact are summarized in Figs. 40 through 49. The 
maximum and residual stresses of these two elements are also 
listed in Table 2. Normal yield stress for the shaft material 
is 35.0x10' Pa (see Section 3.3). The maximum Y-direction 
stresses of elements 841 and 851 are beyond this elastic 
limit, which causes severe plastic deformations, as shown in 
Table 2. As stated previously, the loading from impact of the 
diamond on the rotating shaft was not symmetrical. If the 
shaft (i.e. target) was not rotating, the time history plots 
of stresses for elements 841 and 851 would be the same (such 
as the similarity between Figures 42 and 43) or symmetrical 
with respect to the line of zero stress (Figures 46 and 49). 
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Because the shaft is rotating in the negative Z-direction, 
causing the diamond head (i.e. impactor) to slide toward 
element 851, the time history plots of the stresses for ele- 
ments 841 and 851 are neither the same nor symmetrical with 
respect to the line of zero stress. The irregularity of 
stress time history between time zero and 5x10 seconds in 
Figures 40 through 49 is also caused by the rotation of the 
shaft. 


Table 2: Maximum and Residual Stresses 

for Elements 841 and 851 



Element 

841 

Element 

851 


Maximum 

Residual 

Maximum 

Residual 


( xlO 7 Pa) 

( xlO 7 Pa) 

( xlO 7 Pa) 

( xlO 7 Pa) 

X 

-33.00 

-15.63 

-32.20 

-15.39 

Y 

-46.11 

- 5.28 

-46.67 

- 5.83 

Z 

-31.80 

-15.49 

-32.47 

-15.10 

XY 

+14.92 

- 1.66 

-14.66 

+ 1.86 

YZ 

+14.65 

- 3.65 

+14.88 

- 3.24 

ZX 

- 7.18 

- 2.00 

+ 6.55 

+ 1.82 


Figures 50 through 53 depict the time histories of effec- 
tive plastic strain for elements 831, 841, 851, and 861, re- 
spectively. The plastic deformations of elements 831 and 861 
are much smaller than those of elements 841 and 851. The time 
history plots of plastic strain states for elements other 
than 831, 841, 851, and 861 are all zero (not shown), which 
means that they are still within the elastic range after 
impact. 

The foregoing data were based on a 85 gram applied load. 
Four other loads (35, 50, 185, and 285 grams) were also used 
for the simulation. Table 3 and Fig. 54 show the duration of 
impact versus applied static load. Since the duration of 
response is very short and not easily measured in the labora- 
tory, the data was obtained from results of DYNA3D and not 
compared with experimental results. The heavier the load, the 
longer the duration of contact. 

The scratches produced by impact have a characteristic 
sharp cut as the diamond head scribes the surface. This 



17 


gouging phenomenon pushes material toward one side of the 
scratch developing a ridge. The amount of ridging, which is 
taken to be the difference in height before and after being 
scribed by the diamond head, versus applied static load for 
nodes 100 and 102 (Fig. 35) is plotted in Fig. 55 and listed 
in Table 4. Ridging magnitudes at node 102 are greater than 
those of node 100 by approximately 3%. This is because the 
tangential velocity of the shaft is in the negative Z- 
direction which makes the diamond head pile more material 
toward node 102. 


Table 3: Duration of Impact Versus 

Applied Static Load 


Applied Static Load 
( grams ) 

Duration of Impact 
( microseconds ) 

35 

43 

50 

63 

85 

105 

185 

159 

285 

204 


Table 4: Amount of Ridging at Nodes 100 and 102 

Versus Applied Static Load 


Applied Static Ridging at Ridging at 

Load Node 100 Node 102 

(grams) (micrometers) (micrometers) 


35 

50 

85 

185 

285 


2.0619 

5.9120 

14.8192 

24.0882 

31.1080 


2.1150 

6.2622 

15.2220 

24.8319 

31.9121 
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scratch areas were calculated from computer simulations. 
Their magnitudes versus applied static load were listed in 
Table 5. 


Table 5: Actual and Simulation Scratch Area 

Versus Applied Static Load 


Applied Static 

Actual Scratch 

Simulation Scratch 

Load 

Area 

Area 

( grams ) 

(xlO -4 cm 2 ) 

( xlO -4 cm 2 ) 

35 

2.176 

2.018 

35 

2.209 

— 

50 

— 

2.180 

85 

2.356 

2.297 

85 

2.527 

— 

185 

2.150 

2.398 

285 

2.300 

2.452 


In Figure 56 the triangles represent scratch area data 
from computer simulations and the squares indicate experi- 
mental data. It should be noted that in computer simulation 
the initial conditions (such as initial velocity, point of 
impact, etc.) for each applied load would be exactly the same, 
whereas, in a laboratory environment they were performed manu- 
ally. The maximum percent error between experiment and 
computer simulation of the scratch approximated 10%. 


IV. SUMMARY OF FUTURE OBJECTIVES 


The primary objective for the second year of the research 
program is to change the mesh to simulate a composite mate- 
rial. The mesh developed for the solid material was used to 
verify the results from the finite element analysis with ex- 
perimental data. The next task of the research project is to 
configure a simple cylindrical mesh for a composite material. 
For example, a two ply composite will be modeled first and 
then more complicated structures will be configured by de- 
creasing the shell thickness and increasing the number of 
plies . 
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Once the model has been contructed, it will be fine-tuned 
by comparing it with experimental results. The simulation of 
the multidirectional composite to low velocity impact will 
proceed at a steady rate. The results obtained from this 
stage of the project will be important for understanding the 
unseen damage of delamination in multidirectional composites 
when impacted by a foreign object at low velocities. 
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Figure 4. 80x magnification of a shear crack of an impact 
loaded, Hercules AS4/3501-6 graphite epoxy plate (Ref. 11). 
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Figure 5. Finite element mesh and through-the-thickness stress 
output locations along X-axis. Center of impact point is 
located at the origin and a semicircular spread of the load 
within the contact region is shown (Ref. 12). 
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Figure 6. Shear stress distribution in [90/0/90] lay-up. Refer 
to Fig. 5 for z/t position, where t is the thickness of the 
plate and q Q (GPa) is the applied load (Ref. 12). 
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Figure 10. Finite element mesh of a partially closed crack 
with A as the end point of the closure region S . The hybrid 
element for a closed crack is placed at the cracK tip. Eight- 
node isoparametric elements for anisotropic materials are 
located around the crack-tip element (Ref. 15). 



Figure 11. An open interfacial crack subject to point loads at 
mid-crack surface (Ref. 15). 
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Figure 12. Dynamic stress-intensity factors of Mode I 
delamination (Ref. 21). 
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Figure 13. Dynamic stress-intensity factors of Mode III 
delamination (Ref. 21). 
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Figure 19: Flowchart of Computational Process 
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Figure 20: Flowchart of GMESH Program 
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Front View of One-Half Mesh Geometry 
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Figure 23: Side View of Entire Mesh Geometry 
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Figure 24: 45° View of Entire Mesh Geometry 



Figure 25: Former Symmetrical Mesh Geometry 




Figure 26: Former Mesh Geometry 



Figure 27: Mesh Geometry at 12 Micro-Seconds After 

Initial Impact of Diamond on Rotating Shaft 
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Figure 28: Nesh Geometry at 27 Micro-Seconds After 

Initial Impact of Diamond on Rotating Shaft 
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Figure 29: Mesh Geometry at 117 Micro-Seconds 
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Figure 30: Mesh Geometry at 132 Micro-Seconds at Which 

Rebounding Begins 



Figure 31: Mesh Geometry at 147 Micro-Seconds at Which 

Rebounding Continues 
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Figure 32: Mesh Geometry at 198 Micro-Seconds Showing 

Indented Crater and Bending Of indentor 
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Figure 33: 50° View of Mesh Geometry at 198 Micro-Seconds 
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Figure 34: Top View of Entire Mesh Geometry With Indented 

Crater at 198 Micro-Seconds 
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Fiaure 35: Close-up View of the Impact Region of Interest 
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Figure 37: 


Y-Velocity vs. Time at Nodes 1299 and 101 
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Figure 38: Y-Acceleration vs. Time at Nodes 1299 and 101 
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Figure 40: X-Stress vs. Time for Element 841 
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Figure 41: X-Stress vs. Time for Element 851 
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Figure 43: Y-Stress vs. Time for Element 851 
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Figure 45: Z-Stress vs. Time for Element 851 
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Figure 46: XY-Shear Stress vs. Time for Elements 

841 and 851 
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Figure 47: YZ-Shear Stress vs. Time for Element 841 




Figure 48: YZ-Shear Stress vs. Time for Element 851 
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Figure 49: 
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Figure 54: Duration of Impact vs. Applied Static Load 
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Figure 55 : Amount of Ridging at Nodes 100 and 102 

vs. Applied Static Load 
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Figure 56: Comparison of the scratch area calculated by using 

DYNA3D and measured optically. 
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Appendix A 


Failure Criteria Equations 


Tsai-Wu Criterion [6] 

Fiber failure is assumed to occur when 

F U <T l +2F 12 ff l o t +F 22 02 t +F 66 T2 lt +F l ,T l +F 2 ff t^ 1 (1A) 

where and are normal stresses along and across the fiber direc- 
tion, respectively, and is the shear stress. The following 
strength parameters govern the lamina failure: 
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X' 


22 


-1 

YY' 


(2A) 


66 


1 . 
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F 


-0.5 


12 (XX' YY' ) 1/2 


( 3A) 


where X and X' represent the lamina tensile and compressive 
strengths in the fiber direction, Y and Y' the strengths per- 
pendicular to the fiber direction, and S the lamina shear 
strength. 


Maximum Stress Criteria [7] 

Failure is assumed to occur when any one of the following 
conditions is satisfied: 

a.>X, fiber failure 

c£>Y, transverse tension splitting 

T lt- S ' s ^ear splitting 

and for compressive stresses: 

a.<X', fiber crushing 
Y ' , matrix crushing 
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Appendix B 


Finite-Element Technique for Calculating Elastic Strain Energy 


Release Rates [9] 

The relation between the rate of change of compliance and 
the energy release rate during and incremental crack growth in 
a linear elastic material is given by 

G » 1 p 2 dc (IB) 

1 Ha 


where P is the applied load and C is the compliance of the 
body or the inverse of the stiffness K. 

If C. is the compliance of the body before an incremental 
crack growth and C_ is the compliance after crack extension by 
an amount Aa, then z the rate of change of compliance can be 
written as: 


dC - C, - (2B) 

<Ja A Aa * 

Thus : 

G-lpffl 11 (3B) 

I Aa l R x K 2 J 


where and K 2 are stiffnesses corresponding to and c 2 - 

When the computed energy release rate was equal to the 
critical energy release rate at all node points lying near the 
crack tip, the nodes were split into two and the reaction 
forces were applied to the split nodes. These methods were 
the basis for Mahishi and Adams [ 9 J analysis of delamination 
crack growth in composite laminates. 
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Appendix C 


Laminate Code 

A composite laminate is defined by the following code to 
designate the stacking sequence of ply groups, for example the 
group: 


[0 3 /90 2 /45/_45 3 1 s (1C) 

This code is represented in Figures 57 and 58, and contains 
the following features: 

1. ) Starting from the bottom of the plate, at z»-h/2, the 

first ply group has three plies of 0° orientation; 
followed by the next group with two 90° plies; fol- 
lowed by onle 45° ply; and finally the last group with 
three 45° plies. 

2. ) The subscript s denotes that the laminate is symme- 

tric with respect to the midplane or the z=0 plane. 
The upper half of the laminate is the same as the 
lower half except the stacking sequence is reversed. 

3. ) The laminate is symmetric if 

G(z) = 0 ( — z ) (2C) 

and Qij(z) - Qij(-z) (3C) 

where Qij is the ply material modulus. 
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Figure 57 . Typical stacking sequence of a symmetric laminate 
The laminate code as stated in Equation 1C follows an 
ascending order from the bottom ply (Ref. 6). 



Figure 58. Ply orientations as a function of z. This is 
another representation of Figure 18 (Ref. 6). 


